clear
version 14
set more off
cap log close

mata: mata clear

gl OUT 	"${LOCDRIVE}\Results"
gl TEMP	"${LOCDRIVE}\Tempdata"
gl FIG "${LOCDRIVE}\Figures"

loc r1_hh hhid
loc r1_weight hh_weight
loc r1_indiv sbmemno
loc r2_hh y2_hhid
loc r2_weight y2_weight
loc r2_indiv indidy2
loc r3_hh y3_hhid
loc r3_weight y3_weight
loc r3_indiv indidy3

loc ea ea
loc strata strataid
loc admin0 macroregion
loc admin1 region 		
loc admin2 district  	
loc admin3 ward

loc numeraire other

loc cre 1 /*include correlated random effects*/
loc crelev hh_c /*level of correlated random effects -- normally hh_c or eabase_c or admin3base_c*/ 

loc foodlab "Food"
loc goodlab "Goods"
loc servicelab "Services"
loc otherlab "Other"

* food group labels (for descriptive table)
loc food1lab "SF"
loc food2lab "PNFFV"
loc food3lab "ASF"
loc food4lab "OF"

loc gooddisagglist $gooddisagglist 
loc servicedisagglist $servicedisagglist

forvalues r = 1/3 {
	loc xc_dfl_r`r' = ${XC_DFL_R`r'}
	loc dfl_r`r' = ${DFL_R`r'}
	}
	
log using "${TEMP}\Analysis_Tanzania", replace 



** First make a descriptive table with detailed budget share categories 

loc sumist_food ""
loc sumlist_good ""
loc sumlist_service ""
forvalues f=1/4 {
	loc sumlist_food `sumlist_food' w_food`f'
	}
foreach cat in `gooddisagglist' {
	loc sumlist_good `sumlist_good' wd_`cat'
	}
foreach cat in `servicedisagglist' {
	loc sumlist_service `sumlist_service' wd_`cat'
	}
use "${TEMP}\Tanzania_Panel_estready_all.dta", clear

foreach var in `sumlist_food' `sumlist_good' `sumlist_service' {
	count if `var'==.
	recode `var' (.=0)
	}

eststo clear
bysort round: eststo: estpost summarize w_food `sumlist_food' w_good `sumlist_good' w_service `sumlist_service' [aw=wt_cur]
esttab using "${FIG}\Summary_Tanzania_wdetail.tex", ///
	main(mean) aux(sd) nogaps label nonumbers mtitles(2008-09 2011-12 2013-14) ///
	cells("mean(fmt(3))" "sd(fmt(3))") ///
	tex fragment replace


*********************************************
* The demand estimation borrows heavily from code in
* Tricks with Hicks: The EASI demand system *
* Arthur Lewbel and Krishna Pendakur ********
* 2008, American Economic Review ************
*
* Herein, find Stata code to estimate a demand system with neq equations, nprice prices, 
* ndem demographic characteristics and npowers powers of implicit utility
* This Stata code estimates Lewbel and Pendakur's EASI demand system using approximate 
* OLS estimation and iterated linear 3SLS estimation. Note that iterated linear 3SLS is 
* not formally equivalent to fully nonlinear 3SLS (which does not exist in Stata).  
* However, in some contexts they are asymptotically equivalent (see, e.g., Blundell and
* Robin 1999 and Dominitz and Sherman 2005), and we have verified in our data that 
* coefficients estimated using iterated linear 3SLS are within 0.001 of those 
* estimated using fully nonlinear 3SLS. 
* Code to estimate the fully nonlinear 3SLS/GMM version in TSP is available on request
* from the authors.
* This model includes pz,py,zy interactions. See 'iterated 3sls without pz,py,zy.do' for 
* shorter code to estimate the model without interactions.

use "${TEMP}\Tanzania_Panel_estready_all.dta", clear

** Set number of equations and prices and demographic characteristics and convergence criterion /*#ADJUSTHERE*/
global nf "1" /*number of food groups*/
global nprice "3"
global neq "2"
global neqminus1 "1"
global ndem 1
global npowers "5"
global ncre = 2*$neq
global nshift_core = 17 /*number of demand shifters not including price instruments used as controls which are automatically added if prices considered exogenous*/

global slist ""
global nplist ""
global zlist ""
global zslist ""
global npzlist ""


constraint drop _all


** Set a convergence criterion and choose whether or not to base it on parameters

global conv_crit "0.00000000001"
scalar conv_param=1
scalar conv_y=0


** Data labeling conventions:
* data weights: wgt (replace with 1 if unweighted estimation is desired)
* budget shares: s1 to sneq
* prices: p1 to nprice
* log total expenditures: x
* implicit utility: y, or related names
* demographic characteristics: z1 to zndem

*
loc g1 food
loc g2 good
loc ggood 2
loc g3 service
loc gservice 3 
*

global slist ""
forvalues i=1(1)$nprice {
	g s`i' = w_`g`i''
	global slist $slist s`i'
	}

g p1=lpr_food 
g p1_inst1=lpr_food_i1
g p1_inst2=lpr_food_i2
g p1_inst3=lpr_food_i4
gl p1_endog=1 /*indicates that price `i' is exogenous and instruments will not be used*/

g p`ggood'=lpr_good 
g p`ggood'_inst1=lpr_good_i1
g p`ggood'_inst2=lpr_good_i2
g p`ggood'_inst3=lpr_good_i4
gl p`ggood'_endog=1 /*indicates that price for good is exogenous and instruments will not be used*/ 

g p`gservice'=lpr_service
g p`gservice'_inst1=lpr_service_i1
g p`gservice'_inst2=lpr_service_i2
g p`gservice'_inst3=lpr_service_i4
gl p`gservice'_endog=1 /*indicates that price for service is exogenous and instruments will not be used*/ 

loc nshift_pinst = 0
forvalues i = 1(1)$nprice {
	loc nshift_pinst = `nshift_pinst' + (${p`i'_endog}==0)*3
	}
display "number inst controls as shifters: `nshift_pinst' "

global nshift = ${nshift_core} + `nshift_pinst' /*17 if we do not have any price instruments as controls; add 3 for each endog price*/
display "nshift: ${nshift}"


** set the max matrix size big enough to do constant,y,z,p,zp,yp,yz
global matsize_value=100+$neq*(1+$npowers+$ndem+$neq*(1+$ndem+1)+$nshift)
set matsize $matsize_value


** Normalize the expenditure variable (for estimation)
g x_orig=exp_tot_l
su exp_tot_l
loc mean_log_x = r(mean)
g x=x_orig-`mean_log_x'
su x, detail


** Normalize prices
* generate normalised prices, backup prices (they get deleted), and pAp, pBp
* prices here are all raw (except for the food price which is a fisher index)

global nplist ""
global plist ""

* normalize prices by base period and macroregion
forvalues j=1(1)$nprice {
	su p`j' if round==1 & macroregion==${admin0_base}
	loc pnorm = r(mean)
	replace p`j' = p`j' - `pnorm' 
	la var p`j' "`g`j'' price (normalized)"
	global plist "$plist p`j'"
	}
		
* normalize prices by base good 
su p$nprice if round==1
loc pnorm = r(mean)
forvalues j=1(1)$nprice {
	if `j'!=$nprice {
		g np`j'=p`j'-p$nprice  /*normalize by numeraire good*/
		global nplist "$nplist np`j'"
		}
	}
	
forvalues j=1(1)$neq {
	g np`j'_backup=np`j'
	g Ap`j'=0
	g Bp`j'=0
	}
	
g pAp=0
g pBp=0


** Demographic characterists (demand shifters)

* list demographic characteristics: fill them in, normalize and add them to zlist below
clonevar z1=urban
global zlist "z1"

g zs0=1 /*ones*/
loc zs1 hhsize_adeq
loc zs1_cont = 1 // is continuous

loc zs2 wealthscore 
loc zs2_cont = 1

loc zs3 educ_yrs_head
loc zs3_cont = 1

loc zs4 educ_yrs_max
loc zs4_cont = 1

loc zs5 hh_depshare
loc zs5_cont = 1

loc zs6 pcmh_farm
loc zs6_cont = 0

loc zs7 pcmh_ent
loc zs7_cont = 0

loc zs8 pcmh_market
loc zs8_cont = 0

loc zs9 admin01
loc zs9_cont = 0

loc zs10 admin02
loc zs10_cont = 0

loc zs11 admin04
loc zs11_cont = 0

loc zs12 admin05
loc zs12_cont = 0

loc zs13 admin06
loc zs13_cont = 0

loc zs14 admin07
loc zs14_cont = 0

loc zs15 admin08
loc zs15_cont = 0

loc zs16 r2
loc zs16_cont = 0

loc zs17 r3
loc zs17_cont = 0


loc zscount = ${nshift_core}+1


*Instrument residuals

forvalues i=1(1)$nprice {
	if ${p`i'_endog}==0 {
		cap drop pred1 
		cap drop pred2 
		cap drop pred3
		
		regress p`i'_inst1 p`i' i.round
		predict pred1
		replace p`i'_inst1 = p`i'_inst1 - pred1 
		loc zs`zscount' p`i'_inst1
		loc zs`zscount'_cont = 1
		loc zscount = `zscount'+1
		
		regress p`i'_inst2 p`i' i.round
		predict pred2
		replace p`i'_inst2 = p`i'_inst2 - pred2 		
		loc zs`zscount' p`i'_inst2
		loc zs`zscount'_cont = 1
		loc zscount = `zscount'+1
		
		regress p`i'_inst3 p`i' i.round
		predict pred3
		replace p`i'_inst3 = p`i'_inst3 - pred3 		
		loc zs`zscount' p`i'_inst3
		loc zs`zscount'_cont = 1
		loc zscount = `zscount'+1
		}
	}

forvalues zs=1(1)$nshift  {
    clonevar zs`zs'=`zs`zs''
	global zslist "$zslist zs`zs'"
	}

* standardize continuous demand shifters
forvalues zs = 1(1)$nshift {
	cap su zs`zs', detail
	if _rc==0 {
		if `zs`zs'_cont'==1 {
			loc mean = r(mean)
			loc sd = r(sd)
			replace zs`zs' = (zs`zs' - `mean')/`sd' /*continuous demand shifters are standardized*/
			su zs`zs'
			}
		}
	}

*make pz interactions
global npzlist ""
forvalues j=1(1)$neq {
	forvalues k=1(1)$ndem {
		g np`j'z`k'=np`j'*z`k'	
		global npzlist "$npzlist np`j'z`k'"
	}
}

*make y_stone=x-p'w, and gross instrument, y_tilda=x-p'w^bar

egen mean_x = mean(x)
g y_stone=x-mean_x
g y_tilda=y_stone
forvalues num=1(1)$nprice {
	egen mean_s`num'=mean(s`num') /*, by(round)*/
	replace y_tilda=y_tilda-mean_s`num'*p`num'
	replace y_stone=y_stone-s`num'*p`num'
	}
	
g y=y_stone
g y_inst=y_tilda

su y, detail 
loc yvar_min = r(p1)
loc yvar_max = r(p99)

	
* make list of functions of (implicit) utility, y: fill them in, and add them to ylist below
* alternatively, fill ylist and yinstlist with the appropriate variables and instruments

global ylist ""
global yinstlist ""
global npinstlist "" 
global yzlist ""
global yzinstlist ""
global ynplist ""
global ynpinstlist ""
global creplist ""
global crepinstlist ""
global crepylist ""
global crepyinstlist ""

forvalues j=1(1)$npowers {
	g y`j'=y^`j'
	g y`j'_inst=y_inst^`j'
	global ylist "$ylist y`j'"
	global yinstlist "$yinstlist y`j'_inst"
	}
	
forvalues k=1(1)$ndem {
	g yz`k'=y*z`k'
	g yz`k'_inst=y_inst*z`k'
	global yzlist "$yzlist yz`k'"
	global yzinstlist "$yzinstlist yz`k'_inst"
	}
	
forvalues k=1(1)$neq { 
	g ynp`k'=y*np`k'
	global ynplist "$ynplist ynp`k'"
	}

forvalues k=1(1)$nprice {
	if ${p`k'_endog}==1 { /*interact the price instrument with the y instrument if the price is endogenous*/
	    g ynp`k'_inst1=y_inst*p`k'_inst1
		g ynp`k'_inst2=y_inst*p`k'_inst2
		g ynp`k'_inst3=y_inst*p`k'_inst3
		global ynpinstlist "$ynpinstlist ynp`k'_inst1 ynp`k'_inst2 ynp`k'_inst3"
		global npinstlist "$npinstlist p`k'_inst1 p`k'_inst2 p`k'_inst3"
		}
		else if ${p`k'_endog}!=1 { /*interact price with the y instrument if the price is exogenous*/
			if `k'!=$nprice {
				g ynp`k'_inst1=y_inst*np`k'
				global ynpinstlist "$ynpinstlist ynp`k'_inst1"
				}
			}
	}
	
forvalues g=1/$nprice {
	if `g'!=$nprice {
		bysort `crelev': egen crep`g'=mean(np`g')
		bysort `crelev': egen crepy`g'=mean(ynp`g')
		global creplist "$creplist crep`g'"
		global crepylist "$crepylist crepy`g'"
		}
	if ${p`g'_endog}==1 {
	    bysort `crelev': egen crepy`g'_inst1=mean(ynp`g'_inst1)
	    bysort `crelev': egen crepy`g'_inst2=mean(ynp`g'_inst2)
		bysort `crelev': egen crepy`g'_inst3=mean(ynp`g'_inst3)
		global crepyinstlist "$crepyinstlist crepy`g'_inst1 crepy`g'_inst2 crepy`g'_inst3"
		bysort `crelev': egen crep`g'_inst1=mean(p`g'_inst1)
		bysort `crelev': egen crep`g'_inst2=mean(p`g'_inst2)	
		bysort `crelev': egen crep`g'_inst3=mean(p`g'_inst3)	
		global crepinstlist "$crepinstlist crep`g'_inst1 crep`g'_inst2 crep`g'_inst3"
		}
		else if ${p`g'_endog}!=1 {
			if `g'!=$nprice {
				bysort `crelev': egen crepy`g'_inst1=mean(ynp`g'_inst1)
				global crepyinstlist "$crepyinstlist crepy`g'_inst1"
				}
			}
	}
	

** Set up the equations and put them in a list

global eqlist ""
forvalues num=1(1)$neq {
	global eq`num' "(s`num' $ylist $zlist $yzlist $nplist $ynplist $npzlist zs0 $zslist $creplist $crepylist, noconstant)" 
	macro list eq`num'
	global eqlist "$eqlist \$eq`num'"
	}


** Create linear constraints and put them in a list, called conlist

*create symmetry constraints for price matrix
global conlist ""
forvalues j=1(1)$neq {
	local jplus1=`j'+1
	forvalues k=`jplus1'(1)$neq {
		constraint `j'`k' [s`j']np`k'=[s`k']np`j'	
		global conlist "$conlist `j'`k'"
		}
	}
display "$conlist"
	* 12

*add symmetry constraints for yp interactions
forvalues j=1(1)$neq {
	local jplus1=`j'+1
	forvalues k=`jplus1'(1)$neq {
		constraint `j'`k'0 [s`j']ynp`k'=[s`k']ynp`j'
		global conlist "$conlist `j'`k'0"
		}
	}
display "$conlist"
	* 12 120

*add symmetry constraints for pz interactions
forvalues h=1(1)$ndem {
	forvalues j=1(1)$neq {
		local jplus1=`j'+1
		forvalues k=`jplus1'(1)$neq {
			constraint `j'`k'`h' [s`j']np`k'z`h'=[s`k']np`j'z`h'	
			global conlist "$conlist `j'`k'`h'"
			}
		}
	}
display "$conlist"
	* 12 120 121

g all = 1
g urb = urban==1
g rur = urban==0

loc id_varlist `r1_hh' hh_c eabase_c admin3base_c round wt_pan wt_cur `ea' `strata' `admin1' `admin2' `admin3'

keep `id_varlist' urb rur all $slist x_orig y y_stone y_tilda y_inst $ylist zs0 $zlist $yzlist $plist $nplist $npinstlist $ynplist $npzlist $zslist $creplist $crepylist $yinstlist $ynpinstlist $yzinstlist $crepinstlist $crepyinstlist *_inst1 *_inst2 *_inst3 pAp pBp np*_backup Ap* Bp*

la var urb "Urban"
la var y_stone "Y (log real expenditures)"
la var y_tilda "Y instrument" 

forvalues g=1(1)$nprice {
	loc `g'prop = strproper("`g`g''")
	la var s`g' "``g'prop' budget share"
	la var p`g' "``g'prop' price (log)" 
	}
	
forvalues g=1(1)$neq{
	loc `g'prop = strproper("`g`g''")
	la var np`g' "``g'prop' price (log, normalized)"
	}

tempfile Panel_all_estready_3sls
save `Panel_all_estready_3sls', replace

save "${TEMP}/Panel_all_estready_3sls_subgroups.dta", replace


**********************
* Summarize key vars *
**********************

use "${TEMP}/Panel_all_estready_3sls_subgroups.dta", clear

merge 1:1 hh_c round using "${TEMP}\Tanzania_Panel_estready_all.dta", nogen keep(1 3) keepusing (`zs1' `zs2' `zs3' `zs4' `zs5' `zs6' `zs7' `zs8' )
su $slist $ylist $zlist $yzlist $nplist $ynplist $npzlist `zs1' `zs2' `zs3' `zs4' `zs5' `zs6' `zs7' `zs8' $creplist $crepylist [aw=wt_cur]

eststo clear
bysort round: eststo: estpost summarize y_stone y_tilda $slist $plist $zlist `zs1' `zs2' `zs3' `zs4' `zs5' `zs6' `zs7' `zs8'  [aw=wt_cur]
esttab using "${FIG}\Summary_Tanzania.tex",  ///
	main(mean) aux(sd) nogaps label nonumbers mtitles(2008-09 2011-12 2013-14) tex fragment replace

	
* Make descriptive budget share graphs 
foreach pop in urb rur all {	
	xtile qn_`pop' = y1 if `pop'==1, n(5)

	forvalues q = 1/4 {
		su y1 if qn_`pop'==`q'
		loc q`q'_`pop' = r(max)
		}
		
	drop qn_*
	}
	
forvalues g = 1/$nprice {
	foreach pop in urb rur all {
		lpoly s`g' y1 if `pop'==1 & y1 >= `yvar_min' & y1 <= `yvar_max' [aw=wt_pan], nograph generate(s`g'_`pop'_xgrid s`g'_`pop'_lpoly)
		}
	}
	
forvalues g = 2/$nprice {
	foreach pop in urb rur all {
		loc gm1 = `g'-1
		replace s`g'_`pop'_lpoly = s`g'_`pop'_lpoly + s`gm1'_`pop'_lpoly
		}
	}

foreach pop in urb rur all {
	g s0_`pop'_lpoly = 0
	}
	
	foreach pop in urb rur all {
		graph twoway ///
			(rarea s0_`pop'_lpoly s1_`pop'_lpoly s1_`pop'_xgrid, lwidth(none) fcolor(red%15) ) ///
			(rarea s1_`pop'_lpoly s2_`pop'_lpoly s2_`pop'_xgrid, lwidth(none) fcolor(green%50)) ///
			(rarea s2_`pop'_lpoly s3_`pop'_lpoly s2_`pop'_xgrid, lwidth(none) fcolor(blue%80)) /// 
			||, ///
				xline(`q1_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray)) ///
				xline(`q2_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray)) ///
				xline(`q3_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray)) ///
				xline(`q4_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray))  ///
				yscale(range(0 1)) ylabel(0(.1)1, nogrid)  ///
				title("`pop'") ytitle( "Budget share") xtitle("Log Real Expenditures") ///
				legend(rows(1) order (1 "``g1'lab'" 2 "``g2'lab'" 3 "``g3'lab'")) ///
				name(description_w_`pop', replace) 
		graph export "${FIG}\description_w_`pop'.pdf", replace 
		}
		
drop *lpoly *xgrid

*	


*************************
** Look at first stage **
*************************

use "${TEMP}/Panel_all_estready_3sls_subgroups.dta", clear

loc cv05 = 1.960 /*invttail(`N',0.05/2)*/
loc cv01 = 2.576 /*invttail(`N',0.01/2)*/
loc cv001 = 3.292 /*invttail(`N',0.001/2)*/

forvalues g = 1/$nprice {
	regress p`g' p`g'_inst1 p`g'_inst2 p`g'_inst3 $ylist $zlist $yzlist zs0 zs1-zs${nshift_core}, noconstant
	loc r2_g`g': di %9.3f e(r2_a)
	mat results = r(table)
	test p`g'_inst1=p`g'_inst2=p`g'_inst3=0
	loc fstat_g`g': di %9.3f r(F)
	forvalues i = 1/3 {
		loc p`g'_inst`i'_coef: di %9.3f results[1,`i']
		loc p`g'_inst`i'_coef_se: di %9.3f results[2,`i']
		loc sig_p`g'_inst`i' = (abs(`p`g'_inst`i'_coef'/`p`g'_inst`i'_coef_se')>`cv05') + (abs(`p`g'_inst`i'_coef'/`p`g'_inst`i'_coef_se')>`cv01') + (abs(`p`g'_inst`i'_coef'/`p`g'_inst`i'_coef_se')>`cv001')
		loc stars_p`g'_inst`i' ""
		if `sig_p`g'_inst`i''==1 {
			loc stars_p`g'_inst`i' "*"
			}
			else if `sig_p`g'_inst`i''==2 {
				loc stars_p`g'_inst`i' "**"
				}
				else if `sig_p`g'_inst`i''==3 {
					loc stars_p`g'_inst`i' "***"
					}
		}
	}
	
cap texdoc close 
texdoc init "${FIG}\Prices_firststage_food${p1_endog}good${p`ggood'_endog}service${p`gservice'_endog}", replace force
/*note - if any price is not endogenous then this is not a first stage test because the instruments have been residualized -- ignore this table for food0good0service0*/

tex \hline
tex Dependent & \multicolumn{3}{c}{Coefficients on own price instruments} & F-stat & Adjusted R^2 \\
tex variable: & Price Inst 1 & Price Inst 2 & Price Inst 3 & & \\ \hline
forvalues g=1/$nprice {
	tex ``g`g''lab' price 	& `p`g'_inst1_coef'`stars_p`g'_inst1' & `p`g'_inst2_coef'`stars_p`g'_inst2' & `p`g'_inst3_coef'`stars_p`g'_inst3'  & `fstat_g`g'' & `r2_g`g'' \\
	tex 			& (`p`g'_inst1_coef_se') & (`p`g'_inst2_coef_se') & (`p`g'_inst3_coef_se') & \\ 
	}
tex \hline 
tex Included & \multicolumn{3}{c}{Expenditures, demand shifters,} \\
tex &  \multicolumn{3}{c}{macro-region and survey round fixed effects} \\ \hline

texdoc close


****************
** Estimation **
****************

****** SUREG **********

use "${TEMP}/Panel_all_estready_3sls_subgroups.dta", clear

sureg $eqlist [aweight=wt_cur] ///
		, constr($conlist)  /*Note the "noconstant" is embedded in each equation*/

estimates save "${TEMP}\sureg_all", replace
keep if e(sample)
tempfile sureg_all_data
save `sureg_all_data', all

		
****** 3SLS FIRST STAGE *************

*Model: s = \sum_{r=0}^{r=5} b_r \tilde{y}^r + \mathbf{Czs} + \mathbf{Dz}\tilde{y} + \sum_{l=0}^{L} z_{l}\mathbf{A_{l} p} + \mathbf{Bp}\tilde{y} + \tilde{\varepsilon}
* That is -- y to 5th power, demvars_all, demvars_coreXy, demvars_coreXpvars, pvarsXy
* See p832 of tricks with hicks for constraints:

use "${TEMP}/Panel_all_estready_3sls_subgroups.dta", clear

/*nplist and crepinstlist depending on if price var is endogenous, it will be listed as such with its instruments also listed. if it is exogenous it will not appear in this part of the estimation. If price is exogenous the crep`j'_inst1 variable is just the y inst interacted with the np`j'. If it is endogenous it is also interacted iwth price instrument 1 and there is a crep`j'_inst2.*/
loc pvars_endog ""
loc pvars_exog ""
forvalues j=1(1)$nprice { 
    if ${p`j'_endog}==1 {
		if `j'!=$nprice {
			loc pvars_endog `pvars_endog' np`j' crep`j'
			}
		loc pvars_exog `pvars_exog' p`j'_inst1 p`j'_inst2 p`j'_inst3 crep`j'_inst1 crep`j'_inst2 crep`j'_inst3
		}
	}
	
gl endogvarlist "$ylist $ynplist $yzlist $crepylist `pvars_endog'" 
gl exogvarlist "$yinstlist $ynpinstlist $yzinstlist $crepyinstlist `pvars_exog'"

reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)

estimates save "${TEMP}\reg3approx_all", replace
keep if e(sample)
tempfile reg3_all_data
save `reg3_all_data', all

	
****** 3SLS 2-stage *************

use "${TEMP}/Panel_all_estready_3sls_subgroups.dta", clear

replace y=y_stone
g y_backup=y_stone
g y_old=y_stone
g y_change=0
scalar crit_test=1
scalar iter=0
while crit_test>$conv_crit {
	scalar iter=iter+1
	quietly reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)
	if (iter>1) {		
		matrix params_old=params
		}
	matrix params=e(b)
	quietly replace pAp=0
	quietly replace pBp=0
	quietly replace y_old=y
	quietly replace y_backup=y

	*predict with y=1
	*generate rhs vars,interactions with y=1
	forvalues j=1(1)$npowers {
		quietly replace y`j'=1
		} 
	forvalues j=1(1)$neq {
		quietly replace ynp`j'=np`j'
		tempvar crenpy 	
		bysort `crelev': egen `crenpy' = mean(np`j')
		quietly replace crepy`j'=`crenpy'
		} 
	forvalues j=1(1)$ndem {
		quietly replace yz`j'=z`j'
		} 
		
	*generate predicted values
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y1, equation(s`j')		
		}
		
	*set all p, pz, py to zero
	foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist { 
		quietly replace `yvar'=0
		} 
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y1_p0, equation(s`j')		
		}
	
	*refresh p,pz
	forvalues j=1(1)$neq {
		quietly replace np`j'=np`j'_backup
		tempvar crenp	
		bysort `crelev': egen `crenp' = mean(np`j'_backup)
		quietly replace crep`j' = `crenp'
		forvalues k=1(1)$ndem {
			quietly replace np`j'z`k'=np`j'_backup*z`k'
			}
		}
	
	*generate rhs vars,interactions with y=0
	foreach yvar in $ylist $ynplist $yzlist $crepylist { 
		quietly replace `yvar'=0
		} 
		
	*generate predicted values
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y0, equation(s`j')		
		}
		
	*set all p, pz, py to zero
	foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist { 
		quietly replace `yvar'=0
		} 
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y0_p0, equation(s`j')		
		}

	*refresh p only
	forvalues j=1(1)$neq {
		quietly replace np`j'=np`j'_backup
		}
	
	*fill in pAp and pBp
	forvalues j=1(1)$neq {
		quietly replace Ap`j'=s`j'hat_y0-s`j'hat_y0_p0
		quietly replace pAp=pAp+np`j'*Ap`j'	
		quietly replace Bp`j'=(s`j'hat_y1-s`j'hat_y1_p0)-(s`j'hat_y0-s`j'hat_y0_p0)
		quietly replace pBp=pBp+np`j'*Bp`j'	
		quietly drop s`j'hat_y0 s`j'hat_y0_p0 s`j'hat_y1 s`j'hat_y1_p0
		}

	*round pAp and pBp to the nearest millionth, for easier checking
	quietly replace pAp=int(1000000*pAp+0.5)/1000000
	quietly replace pBp=int(1000000*pBp+0.5)/1000000

	*recalculate y,yz,py,pz
	quietly replace y=(y_stone+0.5*pAp)/(1-0.5*pBp)
	forvalues j=1(1)$npowers {
		quietly replace y`j'=y^`j'
		}
	forvalues j=1(1)$ndem {
		quietly replace yz`j'=y*z`j'
		} 

	*refresh py,pz
	forvalues j=1(1)$neq {
		quietly replace ynp`j'=y*np`j'_backup
		tempvar crenpy 
		bysort `crelev': egen `crenpy' = mean(ynp`j')
		quietly replace crepy`j'=`crenpy'
		tempvar crenp 
		bysort `crelev': egen `crenp' = mean(np`j'_backup)
		quietly replace crep`j'=`crenp'
		forvalues k=1(1)$ndem {
			quietly replace np`j'z`k'=np`j'_backup*z`k'		 
			}
		}

	if (iter>1 & conv_param==1) {		
		matrix params_change=(params-params_old)
		matrix crit_test_mat=(params_change*(params_change'))
		svmat crit_test_mat, names(temp)
		scalar crit_test=temp
		drop temp
		}
	quietly replace y_change=abs(y-y_old)
	quietly summ y_change
	if(conv_y==1) {
		scalar crit_test=r(max)
		}
	display "iteration " iter 
	scalar list crit_test 
	summ y_change y_stone y y_old pAp pBp
	}

	
*now, create the instrument, and its interactions yp and yz
quietly replace y_inst=(y_tilda+0.5*pAp)/(1-0.5*pBp)
forvalues j=1(1)$npowers {
	quietly replace y`j'_inst=y_inst^`j'
	}
forvalues j=1(1)$nprice { 
	if ${p`j'_endog}==1 {
		replace ynp`j'_inst1=y_inst*p`j'_inst1
		replace ynp`j'_inst2=y_inst*p`j'_inst2	
		replace ynp`j'_inst3=y_inst*p`j'_inst3	
		}
		else if ${p`j'_endog}!=1 {
			if `j'!=$nprice{
				replace ynp`j'_inst1 = y_inst*np`j'
				}
			}
	if ${p`j'_endog}==0 {
		if `j'!=$nprice {
			tempvar crenpy	
			bysort `crelev': egen `crenpy' = mean(ynp`j'_inst1)	
			quietly replace crepy`j'_inst1=`crenpy'
			}
		}
	if ${p`j'_endog}==1 {
		tempvar crenpy	
		bysort `crelev': egen `crenpy' = mean(ynp`j'_inst1)	
		quietly replace crepy`j'_inst1=`crenpy'
		tempvar crenpy
		bysort `crelev': egen `crenpy' = mean(ynp`j'_inst2)	
		quietly replace crepy`j'_inst2=`crenpy'
		tempvar crenpy
		bysort `crelev': egen `crenpy' = mean(ynp`j'_inst3)	
		quietly replace crepy`j'_inst3=`crenpy'
		}
	} 
forvalues j=1(1)$ndem {
	replace yz`j'_inst=y_inst*z`j'
	} 
	

*with nice instrument in hand, run three stage least squares on the model, and then iterate to convergence
replace y_old=y
replace y_change=0
scalar iter=0
scalar crit_test=1
while crit_test>$conv_crit {
	scalar iter=iter+1
	quietly reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)

	if (iter>1) {		
		matrix params_old=params
		}
	matrix params=e(b)
	quietly replace pAp=0
	quietly replace pBp=0
	quietly replace y_old=y
	quietly replace y_backup=y

	*predict with y=1
	*generate rhs vars,interactions with y=1
	forvalues j=1(1)$npowers {
		quietly replace y`j'=1
		} 
	forvalues j=1(1)$neq {
		quietly replace ynp`j'=np`j'
		tempvar crenpy 	
		bysort `crelev': egen `crenpy' = mean(np`j')
		quietly replace crepy`j'=`crenpy'
		} 
		
	forvalues j=1(1)$ndem {
		quietly replace yz`j'=z`j'
		} 

	*generate predicted values
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y1, equation(s`j')		
		}

	*set all p, pz, py to zero
	foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist { 
		quietly replace `yvar'=0
		} 

	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y1_p0, equation(s`j')		
		}
	
	*refresh p,pz
	forvalues j=1(1)$neq {
		quietly replace np`j'=np`j'_backup
		tempvar crenp 
		bysort `crelev': egen `crenp' = mean(np`j'_backup)
		quietly replace crep`j' = `crenp'
		forvalues k=1(1)$ndem {
			quietly replace np`j'z`k'=np`j'_backup*z`k'		 
			}
		}
	
	*generate rhs vars,interactions with y=0
	foreach yvar in $ylist $ynplist $yzlist $crepylist { 
		quietly replace `yvar'=0
		} 

		*generate predicted values
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y0, equation(s`j')		
		}

	*set all p, pz, py to zero
	foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist { 
		quietly replace `yvar'=0
		} 
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y0_p0, equation(s`j')		
		}

	*refresh p only
	forvalues j=1(1)$neq {
		quietly replace np`j'=np`j'_backup
		}
	
	*fill in pAp and pBp
	forvalues j=1(1)$neq {
		quietly replace Ap`j'=s`j'hat_y0-s`j'hat_y0_p0
		quietly replace pAp=pAp+np`j'*Ap`j'	
		quietly replace Bp`j'=(s`j'hat_y1-s`j'hat_y1_p0)-(s`j'hat_y0-s`j'hat_y0_p0)
		quietly replace pBp=pBp+np`j'*Bp`j'	
		quietly drop s`j'hat_y0 s`j'hat_y0_p0 s`j'hat_y1 s`j'hat_y1_p0
		}

	*round pAp and pBp to the nearest millionth, for easier checking
	quietly replace pAp=int(1000000*pAp+0.5)/1000000
	quietly replace pBp=int(1000000*pBp+0.5)/1000000

	*recalculate y,yz,py,pz
	quietly replace y=(y_stone+0.5*pAp)/(1-0.5*pBp)
	forvalues j=1(1)$npowers {
		quietly replace y`j'=y^`j'
		}
	forvalues j=1(1)$ndem {
		quietly replace yz`j'=y*z`j'
		} 
		
	*refresh py,pz
	forvalues j=1(1)$neq {
		quietly replace ynp`j'=y*np`j'_backup
		tempvar crenpy 	
		bysort `crelev': egen `crenpy' = mean(ynp`j')
		quietly replace crepy`j'=`crenpy'
		tempvar crenp  
		bysort `crelev': egen `crenp' = mean(np`j'_backup)
		quietly replace crep`j'=`crenp'
		forvalues k=1(1)$ndem {
			quietly replace np`j'z`k'=np`j'_backup*z`k'		 
			}
		}

	if (iter>1 & conv_param==1) {		
		matrix params_change=(params-params_old)
		matrix crit_test_mat=(params_change*(params_change'))
		svmat crit_test_mat, names(temp)
		scalar crit_test=temp
		drop temp
		}
	quietly replace y_change=abs(y-y_old)
	quietly summ y_change
	if(conv_y==1) {
		scalar crit_test=r(max)
		}
	display "iteration " iter 
	scalar list crit_test 
	summ y_change y_stone y y_old pAp pBp
	}

*note that reported standard errors are wrong for iterated estimates
reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)

estimates save "${OUT}\reg3_iterated_all_SL_food${p1_endog}good${p`ggood'_endog}service${p`gservice'_endog}_np${npowers}", replace 
keep if e(sample)


** Save data

cap drop __*

tempfile reg3_iterated_data
save `reg3_iterated_data', all
save "${OUT}\reg3_iterated_data_SL_food${p1_endog}good${p2_endog}service${p`gservice'_endog}_np${npowers}", replace 
